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ABSTRACT 

Molecular hydrogen (H 2 ) is the primary component of the reservoirs of cold, dense gas that fuel 
star formation in our galaxy. While the H2 abundance is ultimately regulated by physical processes 
operating on small scales in the interstellar medium (ISM), observations have revealed a tight corre- 
lation between the ratio of molecular to atomic hydrogen in nearby spiral galaxies and the pressure 
in the mid-plane of their disks. This empirical relation has been used to predict H2 abundances in 
galaxies with potentially very different ISM conditions, such as metal-deficient galaxies at high red- 
shifts. Here, we test the validity of this approach by studying the dependence of the pressure - H 2 
relation on environmental parameters of the ISM. To this end, we follow the formation and destruction 
of H2 explicitly in a suite of hydrodynamical simulations of galaxies with different ISM parameters. 
We find that a pressure - H 2 relation arises naturally in our simulations for a variety of dust-to-gas 
ratios or strengths of the interstellar radiation field in the ISM. Fixing the dust-to-gas ratio and the 
UV radiation field to values measured in the solar neighborhood results in fair agreement with the 
relation observed in nearby galaxies with roughly solar metallicity. However, the parameters (slope 
and normalization) of the pressure - H 2 relation vary in a systematical way with ISM properties. A 
particularly strong trend is the decrease of the normalization of the relation with a lowering of the 
dust-to-gas ratio of the ISM. We show that this trend and other properties of the pressure - H 2 re- 
lation are natural consequences of the transition from atomic to molecular hydrogen with gas surface 
density. 

Subject headings: galaxies: evolution - galaxies: ISM - ISM: molecules - methods: numerical 



1. INTRODUCTION 

The abundance of molecular gas in galaxies is set by 
the complex interplay of various formation and destruc- 
tion processes operating in a highly turbulent medium. 
While many of the individual physical mechanisms are 
relatively well understood, such as the formation of 
molecular hydrogen on dust grains, its photo-dissociation 
by ultraviolet (UV) photons in the Lyman- Werner bands, 
or the i mpor ta nce of dust and H 2 s elf-shi el ding, e.g., 
IDraind (fl978h: IHollenbach fe McKed (fl979T) : ISternbertd 
: iDraine fc BertoldjJT 
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we still lack a realistic 
and coherent picture that links together the formation 
of molecular gas, star formation, the turbulent, multi- 
phase structure of the ISM, and the importance of the 
various feedback channels due to star formation. The 
molecular content of galaxies is a key diagnostic that 
provides insights not only into how galaxies evolve and 
grow their stellar component, but also helps to constrain 
the properties of the physical processes (including feed- 
back) that operate in the ISM. In addition, the modeling 
of the molecular ISM provides a crucial theoretical back- 
ground for the interpretation of molecular gas surveys, 
such as those expected in the near future with the Ata- 
cama Large Millimeter/sub-millimeter Array (ALMA). 
In addition to analytical and numerical mod- 



1 Particle Astrophysics Center, Fermi National Accelerator 
Laboratory, Batavia, IL 60510, USA; feldmann@fnal.gov 

2 Kavli Institute for Cosmological Physics and Enrico Fermi 
Institute, The University of Chicago, Chicago, IL 60637 USA 

3 Illinois Mathematics and Science Academy, 1500 Sullivan 
Rd., Aurora, IL 60506, USA 

4 Department of Astronomy &; Astrophysics, The University 
of Chicago, Chicago, IL 60637 USA 



2009; 



iRobertson fc Kravtsov 
Gnedin et al.l 



2008 



2008 



2010; 



iKrumholz et al.l 

,.. ,20091; IQstriker et"aiT 

I Gnedin fc Kravtsov! 120111 ). empirical correlations 
inferred from observations of nearby galaxies are often 
used to predict the molecular gas abundance. In partic- 
ular, the correlation between the surface density ratio of 
molecular and atomic (HI) hydrog en R mn ] = Sh-,/Sht 
and the mid-plane pressure P BX t (jWong fc BlitzJ l2002t 
iBlitz fc Rosolowskvl 12001 120061? has been included in 
semi-analytical mode ls (e.g.. IDutton fc van den Boschl 
1 20091; iFu et al.l I2010D and numerical simulations (e.g., 
iMurante et al.l I2010D to estimate H 2 mass fractions of 
galaxies and their star formation rates, or in order to 
predict the gl obal evolution of the H 2 baryon content in 
the universe (jObreschkow fc Rawling"sl [20091. Implicit 
in this approach is the assumption that the empirical 
correlation continues to hold for galaxies with ISM 
properties that are potentially different from those 
found in galaxies in the local universe. It is clearly 
crucial to test this as sumption, either observationally 
(Fumagall i et al.l [2010D . or, as is the approach of this 
paper, with the help of numerical simulations that are 
based on a firm theoretical modeling of the microphysics 
of the ISM. 

A further complication is the fact that H 2 masses 
and surface densities are typically not directly acces- 
sible to observations. The kinetic temperature (~ 10 
K) of the bulk of the molecular hydrogen in galaxies is 
too low, and the gas sufficiently shielded from UV ra- 
diation, to populate the excit ed levels of the rotationa l 
ladder at a significant level (jShull fc Beckwithl I1982H . 
Hence, H 2 masses are often inferred from the emission 
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of tracer elements and molecules. In particular, the op- 
tically thick emission from the main isotope of carbon- 
monoxide (CO) serves as a relatively reliable tracer of 
H2 mass, at least under conditions typical of molecular 
cloud s in the Milky W a y (lDickmadll975t iDickman et al.1 
19861 iBloemen et alJ ll986t iMalonev fc Black! 119881: 
Lada et all 11994 iStrong fe Mattoxl Il996t iHunter et al. | 
1997t [Dame et alJ l200lUDraine et all 120071: lA bdo et al.1 
20ld Uckermann et al.ll2012f) ~ However, the conversion 
factor between CO luminosity and H2 mass is expect to 
change systematically with the dust-to-gas ratio a nd with 
the s t rength of the interste l lar radiation field (iWilson 
19951: lArimoto et ail [1991 iBoselli et al.1 l2002t llsraej 



2003 iGkwer fc Mac Low! 120111: iKrumholz et al.1 



2011 



Shettv et alJl201lHLerov et al.ll2011l: iGenzel et al.M 2012: 
Feldmann et all l2012t iNaravanan et all 12012ft . This ef- 
fect needs to be included in the numerical modeling of 
empirical relations that are based on CO observations. 

In this paper we use numerical simulations of galaxies 
in a cosmological framework to study the origin of the 
Pcxt — -Rmoi relation and its dependence on the proper- 
ties of the ISM. Our numerical models compute the local 
HI and H2 abundances in the ISM based on a chemical 
network of well understood formation and destruction 
processes. In addition, we compute the CO emission ex- 
pected from our model galaxies to account for variations 
of the CO — H2 conversion factor. We show that, un- 
der Milky Way like ISM conditions, a P GXt — P mo i rela- 
tion similar to that seen in nearby galaxies arises rather 
naturally. We further show that galaxies with different 
dust-to-gas ratios and/or radiation fields also follow a 
-Pcxt ~ -Rmoi relation, but with changes in the normal- 
ization and the slope. We conclude that H2 abundances 
estimated from a P cxt — P mo i relation are not robust if 
these changes are not taken into account. 

The outline of the paper is as follows. In section [2] we 
describe the set-up of our numerical approach, discuss 
details of the data analysis, and introduce the observa- 
tional data sets that we use to compare with our numer- 
ical models. In the subsequent section [3] we present our 
numerical predictions for the P cxt — R m oi relation and its 
dependence on environmental parameters. We also dis- 
cuss the role of the CO-H2 conversion factor and present 
a physical model that captures many of the properties of 
the P ex t — Pmoi relation. We summarize our results and 
conclude in section 21 

2. METHODOLOGY 

2.1. Simulations 

We use the adapt i ve m esh refinement code ART 
(|Kravtsov et al.l 119971 [2002) to simulate the formation 
and evolution of the baryons and dark matter of a Milky 
Way (MW) sized halo (total mass ~ 10 12 M at z = 0) 
in a 6 Mpc h" 1 cosmological box. We rely on the stan- 
dard "zoom-in" method of embedding the Lagrangian 
region of the halo into layers of lower dark matter res- 
olution to reduce the overall computational cost, while 
capturing the large scale tidal fields correctly ()Katzlll99~il : 
Bertschingcr 2001). The dark matter particle masses are 
mDM = 9 x 10 5 Mq h" 1 in the highest resolution region 
and increase by factors of 8 in subsequent lower resolu- 
tion envelopes. The simulations start from cosmological 
initial conditions with Sl m = 0.3, il\ = 0.7, ilb = 0.043, 



H = 70 km s" 1 Mpc" 1 , and cr 8 = 0.9. 

One of our simulations is run fully self-consistently 
down to z = 1. In contrast, "fixed ISM conditions" (see 
below) are imposed on all the other simulations at z = 4. 
The latter runs are continued for additional 600 Myr be- 
fore they are analyzed. At this time the high-resolution 
Lagrangian region harbors a large disk galaxy sitting in 
a halo of virial mass ~ 4.2 x 10 11 Mq and several lower 
mass galaxies. 

We "fix ISM conditions" (see iGnedin fc Kravtsovl 
120111 ) by imposing a spatially uniform dust-to-gas ratio 
that is independent of the gas metallicity. The gas-to- 
dust ratio is a crucial parameter that affects the dust 
shielding and the formation rates of H2. It also enters 
the CO emission model. We specify the dust-to-gas ratio 
Dmw in units of the dust-to-gas ratio in the solar neigh- 
borhood. Furthermore, we fix the normalization of the 
radiation field at 1000 A to Jo = 10 6 photons cm" 2 s" 1 
sr" 1 eV" 1 , a val ue typical for the solar neighborh ood in 
the Milky Way (IDrainel \WM IMathis et al1ll983h . We 
use the notation Umw — 1, where Umw is the intensity 
of the radiation field at 1000 A in units of Jo- We stress 
that only the normalization of the radiation field is fixed. 
The shape of the radiation spectrum is not modified. 

All simulations include a photo-chemical network that 
follows the formation and destruction of molecular hy- 
drogen in addition to the five major atomic and ionic 
species of hydrogen and helium. The simulations also 
include metal enrichment from supernova (type la and 
type II), but no thermal energy injection, optically thin 
radiative cooling by hydrogen (including H2), helium, 
and metal lines, and 3D radiative transfer of ionizing 
and non-ionizing UV radiation from stellar sources in the 
Optically Thi n Variable Eddington Tensor (OTVET) ap- 
proximation (jGnedin fc Abelll2001h. The details o f the 
implementation can be fou nd in lGnedin et al.1 (|2009l ) and 
IGnedin fc Kravtsovl ([2lnTh . 

We compute the emission arising from the J = 1 —> 
rotational transition of the 12 C 16 Q isotope in a post - 
processing step as described in iFeldmann et all (2012). 
In brief, the CO abundance in each simulation grid 
cell is computed based on the results of a suite of 
small scale magneto-hydr odynamical ISM simulations 
(Gl over fc Mac Low|[2Mll) . The emission is then com- 
puted using the escape probability formalism and assum- 
ing a virial scaling of the CO line width. The contribu- 
tions from these individual, ~ 65 — 130 pc sized resolution 
elements are then combined in the optically thin limit to 
derive the CO emission from larger regions. The CO 
brightness temperature is a free parameter of the model. 
We adopt an excitation temperature of 10 K (approxi- 
mately the kinetic temperature of molecular clouds in the 
Milky Way) and a corresponding brightness temperature 
against the CMB background of 6.65 K. 

An overview of the set of the analyzed simulation snap- 
shots is given in table [TJ 

2.2. Measuring P ext and R mo i 

The ISM pressure in the mid-plane of a disk galaxy is 
not easily accessible to observations. In order to arrive 
at an estimate, disk galaxies are often modeled as self- 
gravitating two-component disks, consisting of a turbu- 
lent gas layer and a stellar disk, in hydrostatic equilib- 
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TABLE 1 

Details of the numerical simulations and the analyzed snapshots 



label 


Ax 


fixed ISM 


%W 


Umw 


rcdshiit 


MW-d-u 


65 pc 


yes 


d = 0.03, 0.1, 0.3, 1,3 


u = 0.1, 1, 10, 100 




MW 


65 pc 


no 






2 = 3 


MW 


130 pc 


no 






z = 1 



Note. — The first two columns provide the name of the simulation and the peak grid resolution of the analyzed snapshots. The 
third column specifies whether dust-to-gas ratio Dmw an d the strength of the UV interstellar radiation field Umw are determined self- 
consistently, or whether "fixed ISM" conditions are used. In the latter case the range of the imposed -Dmw an d Umw values is provided 
in the following two columns. The cosmological redshift of the analyzed snapshots of the self-consistent run is given in the final column. 



rium (jElmegreeni n"989). If these assumptions are made, 
the mid-plane pressure can be expressed in terms of the 
gas and stellar surface densities (T, g , £*) and the velocity 
dispersions of the gas and stellar disk (a g , cr»): 



Pr 



:G£o 



(Tu, 



(1) 



The velocity dispersions can be computed for virial- 
ized disks and equation (fTJ) can be rewritten using the 
scale heights h a and h* of th e gas and stellar disks 
(Blit z fc Rosolowskvll2(M |200(I : 



0.84G°- 5 £ g (Tg 



0.0 



4 h 



0.5 



(2) 



In this paper we follow Blitz & Rosolowski and adopt 
equation @, with Y> g and £* as the independent vari- 
ables. In addition, we adopt their choice a g = 8 km s . 
We also decided to fix the values of the gaseous and stel- 
lar scale heights to h g = 100 pc and h* — 300 pc, typical 
of disk galaxies in the local universe. This is done, be- 
cause in our simulations the disk scale heights are, at 
best, only marginally resolved. However, since the scale 
heights enter in form of a square root, none of our 
results change significantly if the scale heights are varied 
within reasonably bounds. 

H2 masses and surface densities are often derived from 
CO observations using the galactic CO-H 2 conversion 
factor. Its numerical value has been determined to 
within a fac tor of two by a number of independent tech- 
niques fe.g..lSolomon et al.lll987t [Strong fc MattoxllT996l ; 
iDame et all I2001D and a commonly adopted value is 
X CO = N H JWco ~ 2 x 10 20 cm" 2 K" 1 km" 1 s (without 
He), or, equivalently, aco = Mh 2 /Lco ~ 4.4 M pc~ 2 
K _1 km -1 s (including He). 

In order to assess the importance of potential conver- 
sion factor variations, we compute in two different ways 
the H2 surface density that enters the neutral gas surface 



density T, g — 
density ratio 



£h 2 + ^hi (including He) and the surface 



R 



mol 



£h 2 
£hi 



(3) 



The first method takes the H2 surface density directly 
as predicted by the chemical network in the simulation. 
Alternatively, we use the CO emission predicted by our 
CO model (see section I2.1|) and convert it into an H2 
surface density using the galactic conversion factor. In 
general, this has three effects if the actual conversion fac- 
tor differs from the adopted galactic value: (1) it changes 
the minimum H2 surface density that is detectable, (2) 
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Fig. 1. — The P C xt — Rmol relation as predicted by simulations 
that follow the formation and destruction of H2 explicitly. The pan- 
els show the results for a MW-like dust-to-gas ratio (Dmw = 1) 
and for a MW-like UV radiation field (C/mw = !)• The black 
crosses show the pressure inferred from equation (f2j and the ratio 
of the H2 surface mass density (inferred from CO emission with 
a galactic CO-H2 conversion factor) to HI surface mass density 
for 1 kpc 2 patches of the ISM above the CO intensity detection 
limit of 0.2 K km s — 1 . The black solid line is the result of an 
ordinary linear regression of log in ii mo i given log 10 P e xt, in the 
range between the vertical dashed line (the truncation limit) and 
log 10 Pext/(fc B Kcm" 3 ) = 7.5. The fitted slope is 0.69. The trun- 
cation limit is determined as the pressure that corresponds to the 
lowest il mo i value in the sample of kpc 2 patches with a CO emission 
above the CO detection limit. The black dot-dashed line separates 
the regions of data lying below (to the left) or above (to the right) 
the CO detection limit. The red dash ed li nes show the scaling 
predicted by a simple model (see section [3.311 . 

it changes the estimate of P ox t via the change in £ s , and 
(3) it changes i? mo i- 

2.3. Fitting the P cxt — i? mo i relation 

In our set of simulations we measure the P ex t — R mo \ 
relation on kpc 2 patches (the line of sight depth is 1 kpc) 
of the ISM. In order to characterize the main properties 
of the relation, we perform a linear regression in log-log 
space, i.e., we fit for the normalization A and the slope 
S: 



log 10 R mo i = A + S [logj 



5.5]. 



310 fc B Kcm- 3 

We follow iBlitz &: Rosolowskvl (|2006l ) in including in 
our fit only those data points that have P e xt larger than 
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a truncation limit to avoid biasing the regression. In- 
cluding data point at low pressure values can potentially 
bias the slope low, because low P m oi values are preferen- 
tially excluded due to the finite CO detection limit, but 
it can also bias it high, since at i? mo i < 1 the P cxt — i? mo i 
relation steepens due to the relatively rapid transition 
between HI and H 2 , see section [531 

The pressure of the data point that is above the de- 
tection limit and has the lowest P mo i is taken as the 
truncation limit. For instance, in the fixed ISM sim- 
ulation with Umw = 1 an d -Dmw = 1 the trunca- 
tion limit is ~ 1.5 x 10 4 fee K cm" 3 , while in the sim- 
ulation with Umw — 100 and -Dmw = 1 the limit is 
- 9 x 10 4 fc B Kcm- 3 . 

2.4. Observational data 



iBlitz fc Rosolowskvl (|2006[ ) tabulate the results of fits 
to the P ex t — Pmoi relation in their table 2. We take 
the slope and its error bar directly from their table and 
convert the stated normalization Pq to our pivot point 



of log 10 P ox t/(fcB K cnT 



5.5. We exclude the galax- 



ies NGC 598 and NGC 4414 for which we lack reliable 
oxygen abundances. We also exclude their Milky Way 
measurements, because their Fig. 3 shows that most of 
the data falls into the HI-H2 transition regime where a 
si ngle power-l a w is n ot a good fit to the data. 

iLerov et al.l (|2008f) provide in their table 7 radial pro- 
files of the HI, H2 (based on CO observations), and stel- 
lar mass surface densities for a sample of nearby galaxies. 
For each galaxy we compute in each radial bin P ext and 
i? mo i according to ^ and © and perform a linear re- 
gression in the same way as for our simulated galaxies. 
We exclude three of the galaxies (NGC2841, NGC3198, 
and NGC3351) with available H2 data from the fitting, 
either because there are not sufficiently many radial bins 
above the truncation limit to perform a reliable fit or 
because a power-law correlation is not apparent in the 
radial data. 

We note that the binning of the data in radial bins with 
areas >kpc 2 potentially biases the fit results. The reason 
is that the average gas surface density decreases with 
increasing spatial scale due to the inclusion of regions 
that fall below the detection limit at higher resolution. 
Hence, the estimate for P ext decreases and, assuming that 
Pmoi is less affected than P ox t, the normalization of the 
P cx t — Pm ni relation increase s with scale. We therefore 
expect the Lerov et al.l ( 2008) data to be shifted upward 
w.r.t. the Blitz & Rosolowsky (2006) observations. Also, 
since the area of the radial bins varies, the measurement 
of the slope could in principle be affected. However, we 
fi nd that the derived slopes are similar to those found 
in IBlitz fe Rosolowskvl (|2006| ) and, hence, conclude that 
this latter bias cannot be very strong. 

In order to study the P cxt — P m oi relation as function of 
dust-to-gas ratio we supplement the gal axies in both data 
sets w ith gas-phase oxygen abundances (|Moustakas et al.l 
2010) as a proxy for the dust-to-gas ratio. As discussed 
by the authors the absolute calibration of the oxygen 
abundance is rather uncertain (up to ~ 0.6 dex) and de- 
pends crucially on the adopted methodology. We there- 
fore follow the suggestion of the authors and use the 
average of the characteristi c abundances (their table 9) 
derived from a theoret ical dKobulnickv fc Kewlevl 120041) 
and from an empirical (|Pilvugin fc Thuanll2005T) calibra- 



tion method. Half the difference between the two meth- 
ods is used as a measure of the systematic error which 
we then combine with the provided error for each of the 
individual methods to obtain a total error estimate. 

3. THE ORIGIN OF THE i^ X T - Kmol RELATION 
3.1. Predictions of the numerical simulations 
In Fig. Q] we show P ex t, using equation ©, and R mo i, 
using equation ([3]) , for kpc 2 patches of the ISM of a sim- 
ulated galaxy with MW-like ISM conditions. The H 2 
surface density that enters P ex t and P m oi is calculated 
using the CO emission from each patch. The results re- 
main essentially unchanged if we use the actual H2 sur- 
face densities computed in the simulation, because, for a 
MW-like dust-to-gas ratio, the CO-H2 conversion factor 
is close to the canonical galactic value o ver a wide range 
of Hp surface densities (see Fig. 10 in iFeldmann et al.1 
120121 

The figure demonstrates that a power-law relation be- 
tween P GXt and Pmoi arises naturally in simulations that 
follow the microphysics of H2 formation and destruction 
in the ISM. The figure also includes the prediction of a 
simple model which approximates the simulation results 
rather well. We will introduce and discuss this model in 
section 13.31 

In Fig. [5] we compare the slope and normaliza- 
tion of the log 10 P ex t — log 10 Pmoi relation as obtained 
fro m our numerical experime nts w ith the obse r vation s 
by IBlitz fc Rosolowskvl (|2006l ) and ILerov et al.l (f2008l) . 
Overall simulations and observations are in reasonable 
agreement although there are a few interesting outliers. 

Nota bly is the higher no r maliza tion of three galaxies 
in the IBlitz fc Rosolowskvl (|2006f ) sample. As pointed 
out by the authors these three galaxies are in an inter- 
action stage and hence the choice o~ g — 8 km s _1 may 
underestimate the actual gas velocity dispersion. If this 
is indeed so, then P ext is underestimated in these galax- 
ies and, hence, the normalization is overestimated. In- 
terestingly, one of the galaxies in the ILerov et al.l (2008) 
sample, NGC 628 also has a normalization significantly 
above the numerical predictions. However, the values 
of P ex t in the radial bins of this galaxy do not exceed 
log 10 P xt/(fcB K cm~ 3 ) — 5.2 and, hence, the normaliza- 
tion at the pivot point is based on an extrapolation. 

Focusing on the predictions of the "fixed ISM" simula- 
tions, we find that both slope and normalization change 
with dust-to-gas ratio and the UV radiation field in the 
ISM. Consequently, the relation between P ox t and P m oi 
is not universal. In other words, if one wants to predict 
the H2 abundance based on the mid-plane pressure in 
galaxies one needs to take the systematic changes of the 
slope and normalization as function of dust-to-gas ratio 
and UV radiation field into account. The changes are in 
particular 

• a decrease in the normalization with decreasing 
dust-to-gas ratio, 

• an increase in the slope with decreasing dust-to-gas 
ratio, and 

• a decrease in the slope with increasing UV radia- 
tion field. 

We discuss the origin of these trends in section 13.31 
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Fig. 2. — The slope (left panel) and normalization (right panel) of the log 1n P ext — login fimol relatio n as function of metallicity. Black 
circles, x-s, and squares show observational data from Blitz & Rosolowsky (2006) and Leroy ct al. (2008), see legend. The results from the 
"fixed ISM" simulations are shown as lines with error bars for four different values of the UV radiation field (see legend). The magenta 
left-ward pointing triangles show the combined 2 = 3 and z = 1 results of the fully self-consistent run, see table 02 The x-axis of the 
obser vational data refers to gas-phase oxygen abundances derived from nebular emission lines originating in HII regions ( Moustakas ct al. 
2010). In contrast, the x-axis of the simulation data refers to the dust-to-gas ratio which we convert into an oxygen abundance via 



1 

O/H = log 10 



.92. Note this assumes that oxygen abundance and dust-to-gas ratio scale proportional to each other 



12 - 

down to low metallicities. Error bars of the slope and normalization only include fit errors and are est imated via bootstrapping. The 
error bars in the gas-phase oxygen abundance are dominated by systematic uncertainties (see section 12.411 . Slopes and normalizations for 
simulated galaxies with MW-like ISM conditions are in fair agreement with observational estimates for a sample of nearby (many of them 
MW-like) galaxies (see text for a discussion of the outliers). Interestingly, the simulations predict significant changes in the slope and 
normalization with changing ISM parameters. In particular, a decrease in the dust-to-gas ratio is predicted to reduce the normalization of 



the log 10 P c; 



log ln R 



10 -ft-mol 



relation. 



TABLE 2 
Resolution study 



Ax 


method 


slope 


amplitude 


130 pc 


CO 


0.73 ±0.05 


0.77 ±0.02 


65 pc 


CO 


0.69 ±0.02 


0.72 ±0.01 


32 pc 


CO 


0.69 ±0.03 


0.75 ± 0.02 


130 pc 


H 2 


0.71 ±0.03 


0.82 ±0.02 


65 pc 


H 2 


0.67 ±0.02 


0.79 ±0.01 


32 pc 


H 2 


0.69 ±0.03 


0.80 ±0.02 



Note. — Resolution study of a simulation with constrained ISM 
conditions (Umw = 1 Omw = !)• The first column states grid 
size of the highest resolution elements in each simulation. The 
second column denotes whether the H2 surface densities are taken 
directly from the simulation ("H2"), or whether they are estimated 
based on CO intensities calculated with the CO model described 
in section 12.11 The third and fourth column give the slope and 
amplitude of the log in P C xt — login K mo i relation, respectively, and 
also show 1-(T bootstrap errors of the linear fit. The fit results are 
converged for Ax < 130 pc. 



In order to ensure that our predictions are not affected 
by numerical resolution we have re-simulated the run 
with MW-like ISM condition (D MW = I, U M w = 1) 
both at a two times lower and a two times higher spatial 
resolution. Table [2] shows that our numerical predictions 
are converged. 

Fig. [2] also includes the fit results for galaxies taken 
from the z = 1 and z — 3 snapshots of our fully self- 
consistent simulation. Although some differences can be 
seen (e.g., the slightly lower predictions for the normal- 
ization), they follow more or less the predictions for an 
ISM with a MW-like UV radiation field. This is not sur- 



prising since the (volume- weighted) average Umw value 
in these galaxies is typically in the range 0.3-3. Table [3] 
summarizes the relevant global properties of the simu- 
lated galaxies in the self-consistent run. 

3.2. The role of the CO-H2 conversion factor 

It ha s been d emonstrated both observationally, e.g., 
IWilsonl (119951); ILeroy et all (120111). and theo re tically , 
e.g.. iKrumho lz et all (|2011h : iFeldmann et all (|2012l) . 
that the CO-H2 conversion factor increases with decreas- 
ing metallicity or dust-to-gas ratio of a galaxy. Conse- 
quently, the H2 surface mass density derived from CO 
emission will be biased if a constant conversion factor, 
e.g., the canonical galactic value, is used for a galaxy 
with a dust-to-gas ratio or metallicity very different from 
that of the Milky- Way. If so, the estimates for P oxt and 
-Rmoi will be affected and, therefore, the slope and nor- 
malization of the log 10 P C xt — log 10 ^moi relation. In fact, 
the fit results that are shown in Fig. [2] use a galactic 
conversion factor to convert the CO emission into an H2 
surface density. This was done to allow for a fair com- 
pa rison between our simulati ons a nd the obse r vation s 
by iBlitz fc Rosolowskvi ((20061 ) and llerov et all (f2008h . 
However, one may wonder whether, and to which extent, 
some of the trends with dust-to-gas ratio seen in Fig. [2] 
are an artifact of using an incorrect CO-H2 conversion 
factor. 

To address this question we compare in Fig. [3] the 
slopes and normalizations that we obtain using either 
CO-based H 2 surface densities (and assuming a galactic 
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TABLE 3 

Fits to the pressure - H2 relation for simulated galaxies 



Galaxy-ID 


z 


( 1O 1O M ) 


(D M w) 


(^MW) 


slope 


amplitude 


1 


3 


4.2 


0.46 


3.2 


0.86 ±0.31 


0.36 ±0.06 


2 


3 


0.92 


0.29 


0.73 


1.61 ±0.41 


-0.51 ±0.15 


3 


3 


0.33 


0.23 


0.36 


1.67 ±0.82 


-0.83 ±0.23 


4 


3 


0.24 


0.12 


0.30 


1.83 ± 0.43 


-1.11±0.12 


5 


3 


0.22 


0.13 


0.28 


1.40 ±0.57 


-1.06 ±0.17 


6 


1 


11.3 


0.35 


2.4 


1.26 ±0.17 


-0.07 ±0.13 


7 


1 


4.1 


0.56 


1.2 


1.11 ±0.22 


0.13 ±0.11 


8 


1 


0.64 


0.25 


0.15 


0.69 ±0.33 


-0.33 ±0.08 


9 


1 


0.44 


0.23 


0.18 


1.11 ±0.29 


-0.45 ±0.06 


10 


1 


0.28 


0.24 


0.07 


1.05 ±0.31 


-0.33 ±0.07 



Note. — Properties of self-consistently simulated galaxies at two different redshifts. The first column provided the identifier of the 
simulated galaxy. Further column contain: (2) rcdshift, (3) the stellar mass within 10 kpc from the center of the galaxy, (4) the mass- 
weighted average dust-to-gas ratio within 10 kpc (relative to the Milky Way) , (5) the volume-weighted average UV radiation field within 
10 kpc (relative to the Milky Way), (6) the slope of the log 10 P cx t — log 10 R mo i relation, and (7) the amplitude of the relation. The spatial 
resolution is 65 pc at z = 3 and 130 pc at z = 1. 
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Fig. 3. — The impact of the CO-H2 conversion factor on the slope (left panel) and normalization (right panel) of the log 10 P ex t — 
log 10 P mo i relation. Shown are the fit results based on simulations with a fixed dust-to-gas ratio -Dmw converted to oxygen abundance 
via 12 + log 10 O/H = log 10 -Dmw + 8.92 and with a UV radiation field Cmw = 1 (solid lines) or C^mw = 100 (dashed lines). Biases in the 
H2 surface density estimates can affect the fit since both P ex t and P mo i depend on the H2 surface density. Blue lines show the fit results 
if the H2 surface density estimate is based on the CO emission assuming a galactic conversion factor, while red lines use the H2 sur face 
density taken directly from the simulation. The black dot-dashed line shows the prediction of the simple model described in section 13.31 
While some of the trends are aggravated by the dependence of the conversion factor on the dust-to-gas ratio, most of the trends remain 
qualitatively unchanged if the true H2 surface density is used and conversion factor effects are eliminated. 

conversion factor) or by taking the H2 surface densities 
directly from the simulations. The figure shows that the 
decrease in the normalization of the P ex t — Rmo\ relation 
with decreasing metallicity is also present if the true H2 
surface densities are used, showing that it is not primar- 
ily a result of the metallicity dependence of the CO-H2 
conversion factor. The latter does, however, aggravate 
the decrease in normalization at low dust-to-gas ratios. 
The conversion factor is also not responsible for the low- 
ering of the slopes with increased UV radiation field. It 
does, however, contribute to the increase in slope with 
decreasing dust-to-gas ratio. 



3.3. Modeling the P cxt — -Rmol relation 



The results presented in section 13.11 show that a rela- 
tion between P ext and -R mo i arises naturally if H2 forma- 
tion and destruction processes in the ISM are followed 
in a self-consistent manner. We have also demonstrated 
that the numerical predictions are in agreement with ob- 
servational data for nearby, MW-like galaxies and, in 
addition, made specific prediction for the scaling of the 
slope and normalization of the relation with dust-to-gas 
ratio and UV radiation field of a given galaxy. In this sec- 
tion, we will discuss the origin of these trends and show 
that many of them can be understood as a consequence 
of the HI to H2 transition. 

Many observational studies fit the P ex t — -Rmoi relation 
in the H2 dominated regime, i.e., for data with R mo \ > 1. 
In this regime the atomic hydrogen is close to its satu- 



The environmental dependence of the pressure - H2 relation 



7 



ration limit (e.g.. iGnedin fe KravtsovlfeOllD and, hence, 
to first order, P mo i oc E s . Let us now assume that either 
(case A) E*/ft* <C E s //i s (this, for instance, holds in a 
gas-dominated galaxy), or that (case B) E* oc E^. Com- 
bining equations ([2]) and ([3]) we then obtain P m oi oc P^t 
in case A and P mo i oc P^ 2+l3 ' > in case B. 



This simple analysis therefore predicts that P mo i is cor- 
related with P ext and that the slope of the log 10 P e xt — 
log 10 -Rmoi relation in the H2 dominated regime is ~ 0.67 
or 2/(2 + (3) depending on whether the density in the 
mid-plane is dominated by the ISM or the stellar com- 
ponent. In fact, a slope of the order of ~ 0.7 is in- 
deed what we find in our numerical simulations for the 
case Umw ~ 1: see Fig. [T| and EH The obse rvations by 
iBlitz fc Rosolowskvl (|200(f and llerov et all (|2008l ) are 
largely consistent with this slope although with a sub- 
stantial scatter that includes galaxies with significantly 
larger slopes. 

We can quantify these arguments further by using a 
simple model that is bas ed on the numerical s t udy o f the 
HI to H 2 transition by IGnedin fc Kravtsovl (|201lD . It 
is therefore, ultimately, based on the small scale physics 
of H2 formation and destruction in the ISM. It captures 
the basic mechanisms behind the P ex t — Pmoi relation, as 
shown in Fig. [T] and Fig. [3] It works as follows. 

The H2 and HI surface densities are computed as 



(l + E£/(2E,))s 



and Ehi = E„ - Eh 2 



(4) 



where Egj is the HI saturation limit (see below), which is 
a function of dust-to-gas ratio and interstellar UV field. 
With Ehi and Eh 2 given, R mo \ can be computed for any 
given Eg. In order to compute P oxt via equation ([2} an 
estimate of the stellar surface density E* is required. The 
lines shown in Fig.[TJ[3j and IH assum e that Y>* = E„. We 
note that most of the galaxies in the iLerov et"al~l (2008) 
sample show a scaling E* oc E^ with (3 ~ 1 — 2 in the H2 
dominated regime. 
The HI saturation limit is given by the fitting formula 



Eh ! = 20 Mq pc" 



A 4 / 7 



Dmw (1 + Dmw Umw) - 25 ' 



(5) 



In the H2 dominated regime this formula is significantly 
more accurate than the cor responding expressio n (thei r 
equation (14)) provided by IGnedin fc Kravtsovl (|2011f ). 
The A parameter depends on Umw and -Dmw as follows 
(|Gnedin fc KravtsovfeoilD : 



A = In 



(i + . 9 ^(C/mw/15) 4/7 ) , 



a = 



1 + as + s z 



l + s 



0.04 



1.5 x IO- 3 ln(l + (3C/MW) 1 - 7 ) + £>mw ' 

„ . P]\1W 
a = Z.5 ; ; — T77 . 

1 + (C/mw/2) 2 
This model has the following properties: 



The asymptotic slope: 
surface density is 



The ratio of the H? to the HI 



-111 for s s > 



E^j and R 

mol 



Hence, P mo i ~ S 5 /Ej 

4(E g /E H ° I ) 2 for E g < Eg, respectively. Following the 
same argument as outlined in the beginning of this sec- 

tion we infer that P ox t oc Eg . Here, p is 1 if the 
mid-plane density is dominated by the gas component, 
otherwise (3 is the exponent of the assumed E* oc E^ 
scaling. Combining these results we find that in the H2 
dominated regime 



Pmol oc P 



2/(2+/3) 
cxt ! 



i.e. Pmol oc P ° x t 7 for gas dominated galaxies, while in 
the HI dominated regime 



R 



mol 



oc P 



4/(2+0) 



i.e. Pmol oc P cx t f° r gas dominated galaxies. These 
slopes are reasonably close to the slopes ~ 0.8 (H2 dom- 
inated regime) and ~ 1.2 (HI dominated regime) mea- 
sured just a bove and below the transition regime by 
ILerov et all (|2(jol . 

The normalization: — The model reproduces the varia- 
tion of the normalization of the log 10 P ex t — log 10 P m oi 
relation with dust-to-gas ratio as demonstrated by Fig [3l 
This implies that the decrease in the normalization with 
decreasing Dmw is primarily a result of the increase of 
the HI saturation limit with increasing dust-to-gas ratio. 
The model predicts that in the H 2 dominated regime, 
to first order, Ejg oc P^w an d> nen ce, P mo i oc Dmw 
at a fixed E g . We therefore expect the normalization 
of the log 10 P ex t — log 10 Pmoi relation to scale roughly 
linearly with the dust-to-gas ratio, in agreement with 
what our numerical simulations predict. In the HI domi- 
nated regime P m oi oc (E^) -2 and we expected a steeper 
{D MW ) dependence of the normalization on the dust-to- 
gas ratio. 

Observations that convert CO emission into H2 surface 
densities using a constant conversion factor will find an 
even faster decline in the normalization with dust-to-gas 
ratio. A scaling of the conversion factor with dust-to-gas 
ratio of the form P^w implies that the decrease of the 
normalization with dust-to-gas ratio is enhanced by the 
additional factor Pj^w- This effect can be seen quite 
clearly in Fig. [3] In fact, if the model outlined in this 
section is adopted as a baseline, the dust-to-gas ratio de- 
pendence of the conversion factor could be inferred from 
an observational study of how the normalization of the 
log 10 P cxt — log 10 Pmoi relation changes with the dust-to- 
gas ratio (or metallicity) of the galaxy. 

Limitations: — As pointed out, the model predicts the 

scaling oc P cx t oc E g 2+/3 ^ 2 in the H 2 dominated regime 
with /3 ~ 1 if the galaxy is gas dominated. Yet, while 
many of the observations have a slope consistent with 
this prediction, some have a significantly larger slope of 
about unity. In addition, the left panel of Fig. [2] reveals 
that the slope becomes smaller if the UV field is increased 




Fig. 4. — Same as Fig. [T] but for a large UV radiation field (Umw — 100, left panel) and a low dust-to-gas ratio (Dmw = 0.1, right 
panel), respectively. The simple model (see text) fails at reproducing the slopes of the log 10 P cx t — logio Rmol relation recovered via a linear 
regression for the reasons outlined in the text, but it does reproduce the normalization of the data points above the truncation limit. 



and that it becomes larger as the dust-to-gas ratio is 
decreased. How do we explain these trends? 

Formally, a value /? close to zero could explain a slope 
~ 1 in some obs erved galaxies. However, using the 
iLerov et~ahl (|2008|) data we find that (3 lies in the range 
1 — 2 for most galaxies in their sample. We would thus ex- 
pect slopes only in the range ~ 0.5 — 0.7. Instead, the an- 
swer is that in those galaxies the HI surface density starts 
to decrease with increasing S g at large values of the neu- 
tral gas surface density. Note that such a decrease is not 
present in the model, which instead predicts a monotonic 
increase of Shi with S ff , see dU). However, a drop of Shi 
with increasing S g has been observed before, particularly 
in centers of galaxies (|Morris fc Lolll978t iWong fc Blitzl 
12002ft although its origin remains unknown. The decrease 
of Shi with S g boosts the increase of i? mo i with P ext and 
leads to a steeper slope. 

The trend of a decreasing slope with increasing Umw 
appears to be driven largely by the scatter in Shi at a 
given S 9 . Specifically, plotting Shi vs S 9 reveals that 
while Shi does not exceeds its saturation limit, it does 
frequently lie below it. This downward-only scatter can 
be as large as 1 dex at S g ~ 30 M© pc~ 2 and gradually 
decreases with increasing S s out to large surface densi- 
ties (S 9 ~ 1000 M© pc~ 2 ). This change in the scatter 
leads to an overall increase in the median Shi with S g . 
Consequently, i? mo i increases less quickly with P ext and, 
hence, a reduced slope is obtained, see Fig. 

Finally, the increase in the HI saturation limit with 
decreasing dust-to-gas ratio leads to the anti-correlation 
between slope and Dmw- When the dust-to-gas ratio 
is decreased, a higher gas surface density and, hence, 
P cxt needs to be reached in order to fully convert the 
gas from HI to H2. Therefore, a significant fraction of 
the ISM regions with Sh 2 above the detection threshold 
contain a large, or even a dominating, contribution of HI 
and i? mo i > 1 is often not fulfilled. This is illustrated in 
Fig. [4] which shows that cells above the truncation limit 



can have H2 to HI ratios as small as 1%. The figure 
also shows that P mo i has a much steepei0 dependence 
on P ext in the HI dominated regime (i? m oi < 1)- Hence, 
including data points that follow this steeper dependence 
will bias the fit high. A possible solution to avoid this 
problem is to restrict the fit to data that falls either only 
in the H2 or in the HI dominated regime. 

4. SUMMARY AND DISCUSSION 

We studied the relation between mid-plane pressure 
and the Sh 2 /Shi ratio in a set of hydrodynamical simula- 
tions that follow explicitly the formation and destruction 
of molecular hydrogen in the ISM. We have demonstrated 
that these simulation predict a P cxt — Pmoi relation that 
is very similar to the one observed in nearby galaxies. 
The fact that the relation changes systematically with 
the dust-to-gas ratio and the strength of the UV radi- 
ation field in the ISM indicates that the P ex t — P m oi 
relation should not be used as a universal tool to es- 
timate H2 abundances in galaxies with ISM conditions 
different from those in the Milky Way. In particular, we 
find that the normalization of the log 10 -Pcxt — l°g 10 ^mol 
relation decreases with a decreasing dust-to-gas ratio 
Pmw- The scaling of the normalization is cx Dmw &s 
long as the ISM regions that are included in fitting the 
log 10 P cx t — log 10 Pmoi relation are predominantly molec- 
ular and the metallicity dependence of the CO-H2 con- 
version factor is taken into account. We have proposed 
a simple model that is based on the numerical study of 
the HI to H2 transition that captures many of the basic 
properties of the P ex t — P m oi relation. 

5 Empirically it is not too hard to understand what happens at 
Pmoi < !• The surface densities of the star formation rate scales 
super-linear ly wit h S near the so-called threshold of the Kennicut- 
Schmidt relation HSchmidtl fl959l , 119631 ; IKennicuttl 119981) . In this 
regime the atomic hydrogen is not yet fully saturat e d and often 
R mol < 1. Since E S fr oc £h 2 flWone fc Blitz 2002; Bigi el et al.l 
120081) . Sh 2 and, thus, i? mo i scale super-linear ly with S 9 in this 
regime. Hence, i? mo i will also scale super-linearly with P cx t ■ 
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iRobertson fc Kravtsovl ()2008l ) studied the P ext — i? mo i 
relation in hydrodynamical simulations of isolated disk 
galaxies. While their numerical approach differs in de- 
tail, they also find a clear correlation between P e xt and 
i? mo i with a slope ~ 0.9. Although most ISM regions in- 
cluded in their analysis probe the HI dominated regime, 
their largest galaxy reaches pressures high enough to 
study the P ex t — Rmoi relation in molecular hydrogen 
dominated gas. Interestingly, while our numerical mod- 
eling predicts a change in the slope of the relation from 
~ 1.3 for i? mo i < 1 to ~ 0.7 for R mo i J> 1, such 
a change is not immediately apparent in the work by 
IRobertson fc: Kravtsovl (|2008lh As showed in the pre- 
ceding section, a slope of < 0.7 in the H 2 dominated 
regime is a direct consequence of a fixed HI saturation 
lim it. Hence, the relative l y stee p slope (> 0.7) found 
by IRobertson &: Kravtsovl ()2008[ ) in the H 2 dominated 
regime might indicate that in their numerical model Shi 
decreases with increasing E 9 at large gas surface densi- 
ties. 

lElmegreenl (|1993l ) studied the origin of the relation be- 
tween H 2 abundance and mid-plane pressure analytically. 
He found that in the HI dominated regime -Rmoi fn2 
Pext/i) where j is the local radiation field. While this 
relation is sometimes used to "explain" the P oxt — i? mo i 
relation, we stress that the HI dominated regime is not 
the regime typically studied in observations. Further- 
more, j, the interstellar radiation field incident on the 
clouds in the ISM, is likely highly spatially variable and 
will correlate only to a certain degree with galactic prop- 
erties smoothed over ~ kpc scales. This implies that ap- 
plying this formula properly is non-trivial. For instance, 
if the ~ kpc scale UV radiation field is approximately 
constant a naive application of Elmegreen's formula pre- 
dicts i?moi °c -P c 2 x t- In contrast, the model presented in 

section I3~3l predicts R mo \ oc P^} 2+ ^ for i? mo ] < 1, which 
results in a scali ng -R mn i cx P}&__ ~ (for /3 ~ 1 — 2) close 
to observations (jLerov et al.ll2008j ) . 
The numerical predictions presented in this paper are 
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based on the modeling of a variety of well understood 
physical processes in the ISM. However, it is possible 
or even likely that we are missing processes that have 
an impact on the atomic and molecular hydrogen abun- 
dance in galaxies. For instance, while the simulations 
include certain aspects of stellar feedback (e.g., ionizing 
and dissociating radiation from massive stars, metal en- 
richment from supernovae, stellar mass loss), they miss 
others (e.g., radiation pressure, thermal feedback from 
supernovae, stellar winds). Also, the resolution of the 
simulations is too coarse to resolve individual molecular 
clouds and the detailed dynamics and small scale chem- 
istry that takes place within them. Hence, it will be im- 
portant to verify the predictions of our numerical models 
observationally. Important tests include (i) the measure- 
ment of the normalization of the log 10 P ex ± — log 10 -Rmoi 
relation in galaxies with low metallicities and presum- 
ably low dust-to-gas ratios, (ii) the measurement of the 
slope of the relation in both the HI and H2 dominated 
region, and (iii) to check for correlations of the normal- 
ization and slope with the scaling exponent between the 
stellar surface density and the gas surface density. 
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